

/** @file monte_carlo.cc
 *
 * Implementation of monte_carlo 
 *
 */

/*
 *  Copyright (C) 2007 Stefan Weinzierl
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

#include <cmath>
#include <vector>
#include <cstdlib>
#include <iostream>
#include <fstream>
#include <sstream>
#include <ctime>
#include <string>

#include <dlfcn.h>

#include <gsl/gsl_monte.h>
#include <gsl/gsl_monte_vegas.h>

#include "monte_carlo.h"
#include "misc.h"
#include "code_generation.h"

extern "C" double fct_to_integrate(double *u, size_t dim, void *par);

namespace sector_decomposition {

  // -------------------------------------------------
  //
  // Implementation of monte_carlo_parameters
  //

  // ctors

  /**
   *
   * Default constructor
   *
   */
  monte_carlo_parameters::monte_carlo_parameters() :
    iterations_low(), iterations_high(), calls_low(), calls_high()
  {}

  /**
   *
   * Constructor
   *
   * iterations_low is initialised with i_l,
   * iterations_high is initialised with i_h,
   * calls_low is initialised with c_l,
   * calls_high is initialised with c_h.
   *
   */
  monte_carlo_parameters::monte_carlo_parameters(size_t i_l, size_t i_h, size_t c_l, size_t c_h) :
    iterations_low(i_l), iterations_high(i_h), calls_low(c_l), calls_high(c_h)
  {}

  // -------------------------------------------------
  //
  // Implementation of monte_carlo_result
  //

  // ctors

  /**
   *
   * Default constructor
   *
   */
  monte_carlo_result::monte_carlo_result(void) : mean(), error(), chi_squared()
  {}

  /**
   *
   * Standard constructor
   *
   * mean is initialised with mmean,
   * error is initialised with eerror,
   * chi_squared is initialised with cchi_squared.
   *
   */
  monte_carlo_result::monte_carlo_result(double mmean, double eerror, double cchi_squared) : 
    mean(mmean), error(eerror), chi_squared(cchi_squared)
  {}

  // functions 

  /**
   *
   * Returns the mean value.
   *
   */
  double monte_carlo_result::get_mean(void) const
  {
    return mean;
  }

  /**
   *
   * Returns the error estimate.
   *
   */
  double monte_carlo_result::get_error(void) const
  {
    return error;
  }

  /**
   *
   * Returns \f$ \chi^2\f$.
   *
   */
  double monte_carlo_result::get_chi_squared(void) const
  {
    return chi_squared;
  }

  // --------------------------------------------------------
  //
  // Implementation of the Monte Carlo integration
  //

  /**
   *
   * This function is the entry point for the numerical Monte Carlo integration.
   * - expr is the integrand,
   * - x is the list of integration variables,
   * - xmin is the list of the lower integration boundaries,
   * - xmax is the list of upper integration boundaries.
   *
   */
  monte_carlo_result intnum(GiNaC::ex expr, const std::vector<GiNaC::ex> & x, const std::vector<double> & xmin, const std::vector<double> & xmax, const monte_carlo_parameters & mc_parameters, int verbose_level)
  {
    int n = x.size();

    std::vector<GiNaC::symbol> u;
    u.reserve(n);
    std::string name_start("u[");
    std::string name_end("]");

    for (int k=0; k<n; k++)
      {
	u.push_back( GiNaC::symbol( name_start+itos(k)+name_end ) );
      }

    GiNaC::exmap subs_map;
    for (int k=0; k<n; k++) subs_map[ x[k] ] = u[k];
    GiNaC::ex expr2 = expr.subs( subs_map, GiNaC::subs_options::no_pattern );

    double *xl = new double[n];
    double *xu = new double[n];
    for (int k=0; k<n; k++)
      {
	xl[k] = xmin[k];
	xu[k] = xmax[k];
      }

    monte_carlo_result res = intnum_imp(expr2, n, xl, xu, mc_parameters, verbose_level);

    // clean up
    delete[] xl;
    delete[] xu;

    return res;
  }

  /**
   *
   * A simplified entry point for the Monte Carlo integration, where the integration is over the hypercube.
   * - expr is the integrand,
   * - x is the list of integration variables,
   *
   */
  monte_carlo_result intnum(GiNaC::ex expr, const std::vector<GiNaC::ex> & x, const monte_carlo_parameters & mc_parameters, int verbose_level)
  {
    int n = x.size();

    std::vector<double> x_low(n,0.0);
    std::vector<double> x_up(n,1.0);

    return intnum(expr, x, x_low, x_up, mc_parameters, verbose_level);
  }

  /**
   *
   * The implementation of the Monte Carlo integration.
   *
   */
  monte_carlo_result intnum_imp(GiNaC::ex expr, size_t dimension, double *xmin, double *xmax, const monte_carlo_parameters & mc_parameters, int verbose_level)
  {
    typedef double entrypoint(double *uu, size_t ddim, void *ppar);
    entrypoint *module_fct = NULL;

    time_t time_start;
    time_t time_end;
    //if ( verbose_level>0 ) std::cout << std::endl;

    // ----------------------------------------------
    // write the function to a file 
    time(&time_start);
    //if ( verbose_level>0 ) std::cout << " Code generation:       " << std::flush;

    std::ofstream output_file;
    std::ostringstream s;
    srand ( time(NULL) );
    std::ostringstream sss;

    //sss << tempnam("/net/192.168.10.103/hpc_home/users/mkompan/workspace/bogner-sector-decomposition/bogner_tmp/", "xintegrand");

    sss << tempnam("/home/mkompan/workspace/bogner-sector-decomposition/bogner_tmp/", "xintegrand");
    std::string suffix = sss.str();
    s << "" <<  suffix << ".c";
    std::string fName = s.str();
    std::string file_name(fName);

    output_file.open( file_name.c_str() );

    output_file << "// file automatically generated by sector_decomposition " << std::endl;
    output_file << std::endl;
    output_file << "#include <stddef.h> " << std::endl;
    output_file << "#include <stdlib.h> " << std::endl;
    output_file << "#include <math.h> " << std::endl;
    output_file << std::endl;
    output_file << "double fct_to_integrate(double *u, size_t dim, void *par)" << std::endl;
    output_file << "{" << std::endl;

    GiNaC::symbol val_expr("val_expr");
    GiNaC::assign_lst AL(std::string("double"), std::string("t"));
    AL.append(val_expr, expr);
    AL.replace_common_subex();
    AL.split_large_subex(1024);
    output_file << AL << std::endl;
    output_file << " return(val_expr); " << std::endl;

    output_file << "}" << std::endl;

    output_file.close();

    time(&time_end);
    //if ( verbose_level>0 ) std::cout << " Done. Time in seconds: " << difftime(time_end,time_start) << std::endl;

    // ----------------------------------------------
    // compile
    time(&time_start);
    //if ( verbose_level>0 ) std::cout << " Compilation:           " << std::flush;

      std::ostringstream s11;
        s11 << "compile_sector_decomposition " <<  suffix;
        std::string fName11 = s11.str();
    system(fName11.c_str());

    // load the module
    void *module = NULL;
        std::ostringstream s1;
        s1 << "" <<  suffix << ".so";
        std::string fName1 = s1.str();
    module = dlopen(fName1.c_str(), RTLD_NOW);
    if (!module)
      {
	std::cout << "opening of module failed." << std::endl;
      }

    // find the entry point
    module_fct = (entrypoint*) dlsym(module,"fct_to_integrate");

    time(&time_end);
    //if ( verbose_level>0 ) std::cout << " Done. Time in seconds: " << difftime(time_end,time_start) << std::endl;

    // ----------------------------------------------
    // integrate by Monte Carlo
    time(&time_start);
    //if ( verbose_level>0 ) std::cout << " Integration:           " << std::flush;

    double res,err,chi_squared;

    const gsl_rng_type *T;
    gsl_rng *r;
    gsl_rng_env_setup();
    T = gsl_rng_default;
    r = gsl_rng_alloc(T);

    gsl_monte_function F = { module_fct, dimension, NULL};

    gsl_monte_vegas_state *mc_state = gsl_monte_vegas_alloc(dimension);

    // warm up to fix the grid
    mc_state->stage = 0;
    mc_state->iterations = mc_parameters.iterations_low; 
    gsl_monte_vegas_integrate(&F, xmin, xmax, dimension, 
			      mc_parameters.calls_low, r, mc_state, 
			      &res, &err);

    // the main run
    mc_state->stage = 1;
    mc_state->iterations = mc_parameters.iterations_high; 
    gsl_monte_vegas_integrate(&F, xmin, xmax, dimension, 
			      mc_parameters.calls_high, r, mc_state, 
			      &res, &err);
    chi_squared = mc_state->chisq;

    time(&time_end);
    //if ( verbose_level>0 ) std::cout << " Done. Time in seconds: " << difftime(time_end,time_start) << std::endl;

    // ----------------------------------------------
    // clean up
    gsl_monte_vegas_free(mc_state);

    dlclose(module);

       std::ostringstream s2;
        s2 << "rm " <<  suffix << ".c";
        std::string fName2 = s2.str();
    system(fName2.c_str());
           std::ostringstream s3;
            s3 << "rm " <<  suffix << ".o";
            std::string fName3 = s3.str();
    system(fName3.c_str());
           std::ostringstream s4;
            s4 << "rm " <<  suffix << ".so";
            std::string fName4 = s4.str();
    system(fName4.c_str());

    //if ( verbose_level>0 ) std::cout << "                        " << std::flush;

    return monte_carlo_result(res, err, chi_squared);
  }


} // namespace sector_decomposition

